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Violent relaxation during the collapse of a galaxy halo is known to be incomplete 
in realistic cases such as cosmological infall or mergers. We adopt a physical picture 
of strong but short lived interactions between potential fluctuations and particle 



orbits, using the broad framework outlined by Tremaine (1987) for incorporating 
incompleteness of the relaxation. We are guided by results from plasma physics, viz. 
the quasilinear theory of Landau damping, but allow for significant differences in our 



case. Crucially, wave particle scattering does not drive the system to an equilibrium 
distribution function of the exponential type, even in regions of phase space allowed 
by the constraints. The physical process is mixing without friction in "action" space, 
for which the simplest possible model is a constant phase space density modulated by 
the constraints. Our distribution function does not use the exponential functions of 
the energy prevalent in other work, which we regard as inappropriate to collisionless 
systems. The dynamical constraint of a finite short period of the relaxation process 
is argued to lead to a 1/T r factor in the distribution function, where T r is the radial 
period. The notion of strong potential fluctuations in a core is built in as a cutoff in 
pericenter (which we find preferable to one angular momentum, the other alternative 
we explored). The halo of the self-consistent, parameter-free solutions show an r~ 4 
behavior in density at large r, an r 1 ' 4 surface brightness profile in the region 0.1 — 8 r e , 
and a radially anisotropic velocity dispersion profile outside an isotropic core. The 
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energy distribution seen in simulations N(E) singles out the pericenter cutoff model 
as the most realistic among the variants we have explored. The results are robust to 
modifications of the period dependence keeping the same asymptotic behaviour, or to 
the use of binding energy raised to the power of three halves in place of 1/T r . 



Subject headings: galaxies: formation — galaxies: structure — galaxies: interactions, 
violent relaxation 



1. A model of quasi- violent relaxation 

Elliptical galaxies are expected to have undergone violent relaxation (Lynden-Bell 1967), 
which is a collisionless process whereby the energies and angular momenta of stars and dark 
matter particles get redistributed by strong potential fluctuations in such a way that the outcome 
depends mainly on the macroscopic features of the initial conditions. This concept is supported 
by N-Body experiments where for a range of initial conditions, the final state has more or less a 
universal profile. 

One of the insights gained from cold collapse simulations is that a compact region or a core 
develops, the system partially reexpands and then in a few crossing times settles into a centrally 
peaked configuration with most of the mass inside a radius within a fraction of the initial size 
(van Albada 1982). The potential fluctuations are initially of large amplitude but they damp out 
in a few crossing times. Clumpy and cold initial conditions with small values of 2T/W ~ 0.1, 
where T and W are total kinetic and potential energies, are preferred over hot or homogeneous 
conditions in order to produce final configurations that fit the de Vaucoleurs exp(— r 1 ' 4 ) surface 
density profile better. 

The general violent relaxation problem is difficult because the fluctuations have large 
amplitudes and damp in a few crossing times. Thus there is no small parameter and just a 
single time scale. The physical process belongs to the broad category of interactions between 
waves (potential fluctuations) and particles (stars or dark matter), wherein the amplitude of 
the wave is self-consistently determined by the positions of the particles. We start with the 
limiting case of small amplitudes, and long-lived waves, where we might expect the problem to 
submit to perturbative methods. However, it is difficult to determine even the linear modes 
of self-gravitating, collisionless systems; this is because galaxies are spatially inhomogeneous 
objects. Quasi-neutral plasmas are not subject to this problem, so let us examine the quasi-linear 
theory of Landau damping of electrostatic waves in collisionless plasmas (c.f. § 49 of Lifshitz 
& Pitaevskii 1981). The focus of the theory is on the slow evolution of the coarse-grained DF 
(distribution function) of electrons, f(p,t), where p is momentum. A key result is that f(p,t) 
obeys a diffusion equation in p-space: 



m-Wa [ Da "df,J ' (1) 

where the diffusion coefficient, D a g, is nonzero only in a range of p corresponding to electrons 
that are resonant with the perturbations. This diffusion causes smoothing, and creates a plateau 
in f(p). Furthermore, D a p is proportional to the energy in the fluctuations; as this energy is 
absorbed by the resonant electrons, the diffusion itself weakens, and f(p) reaches a steady state. 
The point we wish to make here is that this evidently non-Maxwellian steady state has been 
reached by a self-regulating diffusive process. 

There is of course a basic difference between slow damping of plasma waves, and the damping 
of the oscillations of a violently relaxing galaxy. In the former, the energy in the waves is absorbed 
by resonant electrons, whereas the very concept of resonance is dubious in the latter case, because 
the oscillations are short-lived. In a relaxing galaxy, the orbits of particles are expected to be 
scattered by an oscillating core (Tremaine 1987, Spergel & Hernquist 1992 (SH)). In realistic 
situations, the deflections are individually large but small in number. Let us, nevertheless, for 
orientation, initially consider the limit in which each star, or dark matter particle, undergoes many 
small kicks. Hence particles move under the combined actions of integrable forces, and small kicks 
that lead to global stochasticity. In stochastic regions of phase space far from islands, the diffusion 
of actions is well described by a Fokker-Planck (FP) equation. A feature of the FP equation for 
Hamiltonian kicks is the complete absence of dynamical friction, a result as old as Landau 1937 
(c. f. § 5.4 of Lichtenberg & Lieberman 1983). The FP equation resembles equation (|l|), where p 
are the actions. Once again, D a p is proportional to the square of the perturbation. The absence 
of friction implies that a Maxwellian distribution is not the final steady state. In fact, the process 
being purely diffusive, it is evident that the DF would approach one that is independent of the 
actions, given a finite volume of phase space and enough time. 

We learn from the above examples that in the limit of long-lived, small amplitude fluctuations, 
the relaxation of a galaxy is likely to be primarily diffusive in phase space. This relaxation process 
is very different from collisional relaxation in neutral gases, which is a reflection of the long-range 
nature of gravitational forces in contrast to the short range of interatomic forces. The distinction 
is made clearer by using the H-functions introduced by Tremaine, Henon, and Lynden-Bell (1986). 
They are actually functionals of the coarse-grained DF, defined as 

H[f} = -jc(f)d 3 xd 3 v, (2) 

where C(f) is a convex function. Any mixing process conserving fine grained phase density, 
such as collisionless relaxation, will increase H[f}. This result needs the assumption that the 
initial state had the fine and coarse grained distribution equal to each other, which applies to the 
cold collapse simulations with which we are concerned with. It does not imply that H increases 
monotonically (Dejonghe 1987, Sridhar 1987). On the other hand, binary collisional processes, 



such as those relevant to thermal relaxation in neutral gases, do not conserve phase density and 
single out a unique H-function, namely the Boltzmann entropy given by C{f) = /In/. 

We now extrapolate to the case of the large amplitude, short-lived fluctuations appropriate to 
violent relaxation. When the "actions" change by large amounts, the mixing process is no longer 
correctly described by a diffusion equation. However, the relaxation is probably well described by 
a process of mixing without friction, and an initial DF will spread as far as it can in phase space; 
the extent and nature of the spreading being determined by dynamical constraints discussed in the 
following section. We write the DF = A(£) x factor expressing constraints, where £ is the single 
particle binding energy (bound particles have positive £). In the short duration before freeze-out, 
A spreads out nearly evenly and beyond £ = (some particles escape). But in constructing 
distribution functions of the galaxy we do not include unbound particles. The parallel to diffusion 
without friction suggests a final state for which A(£) = constant . 

Note that this hypothesis is the simplest choice consistent with the physics of violent 
relaxation. A similar concept was proposed by Yan'kov (1994) to explain turbulent transport 
in plasmas of tokomaks with a uniform distribution of particles in a region specified by certain 
constraints of the problem. 



2. A prescription for f(£, J 2 ) 

The relevant dynamical constraints will clearly depend on the context; for instance, the merger 
of two galaxies is very likely to be different from the formation of dark halos by cosmological infall 
over extended periods of time. We now seek the dynamical constraints appropriate to situations 
such as the cold collapses familiar from numerical simulations (c.f. van Albada 1982) which might 
also carry over to the case of head-on mergers between galaxies. For simplicity, we assume that the 
relaxed system is spherical, and its DF is a function of the energy, £, and the angular momentum, 
J. Tremaine (1987) has suggested the following two plausible physical requirements (see also 
Merritt, Tremaine, & Johnstone 1989 (MTJ), Stiavelli & Bertin 1987). 

Rl. The potential fluctuations last for only a limited time, T e , which is of order a few crossing 
times. Hence orbits with radial periods, T r , exceeding T e will be underpopulated by a factor 
~ (T e /T r ). In terms of a plausible physical requirement, this is viewed as a uniform filling of those 
orbital phases which allowed the particle to visit the core in the time T e , and zero elsewhere. 
Finally, of course, all phases are equally populated in the coarse grained function, which means a 
filling proportional to T e /T r . 

R2. The fluctuations are largely confined to a central region or a core of radius r c . Hence only 
orbits whose pericenters are smaller than ~ r c would have visited the region of large potential 
fluctuations, and have had a chance to relax violently. The validity of this will depend strongly 
on the initial conditions; the collapse of a cold, nearly non-rotating initial configuration is more 
relevant to this study. 



With this motivation, we assume the following form for /(£, J 2 ) 

(A function proportional \ / Cutoff in pericenter 
to 1/T r x or | : A (■>,) 

for large T r J \ in J 2 

The first two factors in the distribution function are the dynamical constraints given by Rl 
and R2 respectively, and the third factor, A, is taken to be a constant as discussed earlier. The 
first factor 1/T r ensures the continuity of the DF at £ = 0, i.e., f{£ , J 2 ) = for £ < 0. Jaffe 
(1987) remarked that one asymptotic property of violently relaxed systems is that N{£) should 
be finite and non-zero at £ = because the number of particles ejected from the core should be 
smooth across £ = 0. We show in §Q that for small £, the restricted density of states due to the 
second factor in (||) is ~ £~ 3 ' 2 , while / ~ 1/T r ~ £ 3 ' 2 (for finite mass systems) and therefore 
N(£)~£°. We define 

f(£,J 2 )^f (£,J 2 )-C-A (4) 

where /o and C represent the first and second factors in (g). In this work, we studied the following 
closely related models based on (||). 

• A model with fo(£) = £ i ' 2 was first obtained and used as an input to generate the 1/T r 
counterpart discussed below. Details of this model are presented in appendix |A]. The 
dependence of the distribution function on £ and J is explicit, and hence the analytic work 
can be carried further. 

• fo(£,J 2 ) = 1/T r , is presented in §||. Here, the dependence of the radial period on £ and 
J has to be derived iteratively from the potential. The models above form a good starting 
point for successive numerical iterations which we refer to as r c and J m models. A variant of 
these models with f (£, J 2 ) = l/(T r + T ), with an extra parameter T introduced to assess 
the sensitivity of the results to the functional form of the period cutoff. 

All the above models have the asymptotic properties p ~ r~ 4 , N(£) ~ £° (described above), 
and surface densities approximating the r 1 ' 4 law. The form of energy distributions, N(£), at 
distances larger than the core radius r c , or the scale radius rj, are derived for the pericenter model 



in §|j. In §5.1, we describe the form of N(£) for the T r models. This is helpful in picking out the 
pericenter cutoff models as a better fit to results of simulations. In §|3| we describe the boundaries 
in the £ — J 2 plane required in doing the integrals in the allowed velocity space. The properties of 
the solutions and comparison to simulations are described in §||, and we present a discussion in §|7], 
and conclusions in SS The reader interested in first examining the results is urged to go to 
and then return to §0-0 for some details of the construction. 



3. Pericenter and angular momentum constraints in £ — J 2 plane 

For a spherical model, with a distribution function, f(£,J 2 ) = Afo(£, J 2 ) C, including a 
sharp pericenter cutoff at r c or in angular momentum, J m , Poisson's equation can be written as 

^§) = ^ 2 Jd£ dJ 2 (2(-£ - *) - J>T 1/2 US, J 2 ) C( V ), (5) 

where A is a constant that is determined by normalization. The cutoff function, C, is given by 



Cfa) = ' " (6) 




where rj equals either r p /r c where r p (£, J 2 ) is the pericenter, and r c is the cutoff radius, or J/J m , 
where J m is the maximum angular momentum. We explore both these constraints here and 
construct a method to solve the self-consistent models with f (£,J 2 ) = 1/T r ; eqn (H), coupled 
with 

T r (£,J 2 ) = 2 f ra dr{2{-£-^)-J 2 /r 2 )- 1 / 2 , (7) 



where r p and r a are the turning points. We have solved the coupled equations by a semi-analytical 
method whose details are presented in §|5[ 

Now we introduce a variable, ro(J 2 ), which is the radius of a circular orbit for a given angular 
momentum, J 2 = Tq&q. This is useful in determining the region of integration in the £ — J 2 plane. 
Note that J 2 = M (ro)ro from Poisson's equation, where M(r) is the mass inside r and hence it is 
a monotonically increasing function of r$. Fig. |l] shows a plot of the absolute value of the effective 
potential, — £f(r, ro) = $ + J 2 /(2r 2 ) = $ + &' ro 3 /(2r 2 ) for a given J 2 and a typical <1>. We now 
consider the region of integration allowed by C(r/). 



3.1. Case of 7] = r p /r c 

The region of integration in the £ — J 2 plane is bounded by following curves which are shown 



m 



Fig. § 



Kl. The upper bound of the region of interest is given by the minimum of the effective potential 
(—£f) of a circular orbit for a given $(r) and ro ie, £ < £f(ro,ro) = — ^o — 3>o r o/2 fo r r p < tq < r c . 
The orbits with r p < r c < ro, obey the bound, £ < £f(r c ,ro) = —<& c — < ^ ) o r o 3 /(2r^) as indicated in 
Fig0. 

K2. At a given radius r for the system, the effective potential is bounded by 
£ < B(r, J 2 ) = £j(r,r ) = -$ - J 2 /(2r 2 ) or equivalently v 2 (r) > 0. 

For r < r c , the pericenter constraint is satisfied and the operative bound is just £ < £f{r,ro). 
Hence, the bounding line £ = £f(r c ,ro) given by constraint K2 lies inside and is tangent to the 
curve given by constraint Kl; the point of contact represents a circular orbit for a given ro- 



Now consider the case r > r c . The point of intersection, (£*, J 2 ), between the bounding line of 
bound K2, £ = £f(r, tq) and £ = £f(r c ,ro) given by bound Kl represents an orbit whose turning 
points are r and r c (see Fig ^). The point of intersection works out to be 

£* = ($ c r c 2 - $>r 2 )/{r 2 - r 2 ) (8) 

J 2 = 2(<l>-cl> c )r 2 r 2 /(r 2 -r 2 ) (9) 

Effectively, for r > r c , the bound given by constraint K2, £ < £f(r,ro), applies upto 
J 2 = r^' c , beyond which the bound, £ < — <1> C — J 2 /(2r 2 ) given by Kl is operative, and forms a 
wedge shaped region. For r < r c , the region of integration is given by the bound £ < £f(r, tq) and 
is a triangular shaped region. The regions of integration are shown in Fig. y and are summarized 
by 



At = £ < -$> - J 2 /(2r 2 ), r<r c (10) 

£ < -$ - J 2 /(2r 2 ), r>r c k J 2 < J 2 
£ < -$ c - J 2 /(2r c 2 ) r>r c k J 2 > J 2 

Now that the regions of integration have been determined, we can write the integral on the RHS 
of (pi), without the normalization constant, A as 



- 42=< ! c. * t2/,o„.2n „^„ , r 72^ 7 2 ( U ) 



T , n j 2To(r /o) r < r c 

[ 2b(r;/o)-2L(r;/o)+X + (r;/ ) r > r c 

where J 2 (r) = — 2<l>r 2 is the intercept on the J 2 axis, and 

X (r; /o) = \ \' h dJ 2 / B ° '" ] d£ (2(-£ - $) - J 2 / r 2 )' 1 / 2 f {£ , J 2 ) (13) 

T_(r; / ) = 1 /^ dJ 2 / S(r ' J } d£ (2(-£ - $) - J 2 /V 2 )- 1/2 /o(£, J 2 ) (14) 

r^ 7j2 Jo 

^ + (r; /o) = - / d£ dJ 2 (2(-£ - $) - J 2 /r 2 )- l ' 2 U£, J 2 ) (15) 

H Jo JJl 

3.2. Case of r/ = J/J m 



The constraint K2 applies here and the maximum allowed angular momentum is J m . There 
are no circular orbits allowed outside the radius t given by J^ = — 2<F ; r 2 = £ 3 <J>'(£). Also, for 
r < rj only constraint K2 applies and for r > rj, the orbits with J 2 > J^ are excluded. Hence the 
region of integration, B is given by 



Bi = £ < -$ - J 2 /(2r 2 ) r < r^ 

£ 2 = £<-&- J 2 /(2r 2 ) kJ 2 <J 2 m r> Tj 



(16) 



For r\ = J/J m , the integral on the RHS of (||), without the normalization constant A will reduce to 

Mr;/o) = T , n r r ^ ^ ( 17 ) 

[ M\r\h) -/C-(r;/o) r > r 3 - 

/C_(r;/ ) = 4 /^dJ 2 /*""' ) d£:2(- ( S-cI>)-J 2 /r 2 )- 1 / 2 /o(«?,J 2 ). (18) 



where 



4. Restricted density of states for a pericenter cutoff 

Here we calculate the density of states for a model that has a sharp pericenter cutoff. The 
restricted density of states is given by 

g{£) = J d 3 r C(r p /r c )^ = 8vr 2 J dr j'™ dJ 2 /^2(-£ - $) - J 2 /r^ 

= 16ir 2 y/2J rdr f^ {-£ - <S>)r* - y/(-£ - <^)r^ - J^J2J (19) 

where Jm ax is determined by constraints. Hereafter we work with units where GM = 1. There is 
no contribution to g{£) from £ > — $ and therefore for a given energy, £ , only the region r < rg is 
accessible where $>(rg) = —£. Consult Fig. ^ for the allowed range of integration, A. For r < r c , 
the pericenter constraint (Kl) is satisfied and therefore, J^ax = ^e = 2(— £ — &)r 2 as given by 
K2. Now, take the case when £ > £ c and r > r c shown in Fig. |5| where £ c = —Q c — &' c r c /2, the 
energy of the circular orbit of radius r c . Earlier, the point of intersection of — <£ — J 2 /(2r 2 ) and 
— Q c — J 2 /(2r 2 ) was defined to be £*, J 2 . Here, £ > £ c > £*, and hence the pericenters lie inside r c 
and Jm ax = J e as given by constraint K2. Now consider the case £ < £ c and r > r c as illustrated 
in Fig. | for which J 2 = -2(5 - & c )r 2 < J*- It is clear that J^ ax = Min(J 2 , J 2 ) and if J 2 < J 2 , 
the particle has a Jm ax = J 2 because any higher angular momentum will include orbits whose 
pericenters lie outside r c . This happens only for r c < r < rj where rj(£,r c ) is the apocenter of an 
orbit for a given energy —£ and pericenter r c . It is given by 



r I~ r c 



which has two roots; we seek the one for which 77 > r c . To summarize: 

{(-£ - $)r 2 ; -3> > £ > £ c , r < rg 
(-£ - $ c )r 2 ; £ c >£*>£, r c <r<rj<r £ (21) 

(-£ - $)r 2 ; £ c > £ > £*, r c < 77 < r < r £ 

We may then write 

g{£) = 167r 2 V2 (-g c {£) + f r 2 y/-£-$ dr) (22) 



where 

(FS _ / Ir 1 dr r 2 j£7=£(A ~ S) 1/2 £ c > £* > £ r c <r<r I <r e 
^ U otherwise 

We can now use the equation above to calculate the density of states and the energy distribution 

for a Keplerian potential, the asymptotics of which are applicable more generally for finite mass 

systems. After some algebra, the restricted Keplerian density of states (in units where GM = 1) 

is given by 

f £-5/2 £ > £ 

*< £) = ^{ -4 £-3/^(1-^1) s<S c (24> 

where it is continuous at £ c = l/(2r c ). The energy distribution for fo(£) = A'£ s > 2 , will then go as 



£- 1 £>£ c 

4$ c r 2 (l-£$ c r 2 ) £ < £ c 



N k (£)=A>^V2{ ^ 2n ^ a . t", (25) 



Note that for any finite mass potential, the above equations are valid for large radius 
(r S> r c ) and only approximate in the inner region where the realistic potential deviates from 
1/r. Clearly, as r c — ► oo, £ c — ► 0, and one recovers the Keplerian form of g{£) oc £~ 5 ' 2 and if r c 
is zero then g{£) vanishes as expected. Near small £, the velocities become more radial and the 
unconstrained Keplerian density of states, £~ 5 ' 2 , is reduced by a factor v 2 oc £ . In our model the 
assumption of a sharp pericenter cutoff at r c leads to g{£ ) oc £~ 3 ' 2 , near £ = and the choice of 
f(£) = 1/T r ~ £ 3 ' 2 , at small £ , based on the dynamical arguments made earlier, is consistent with 
the required property of the "break of N(£ = 0)", or the finite and non-zero value of N{£ = 0). 
Jaffe(1987) made the interesting point that the demand of a break in N(£) = f(£)g(£) at £ = 
for a (unrestricted) Keplerian density of states, leads to f(£) ~ £ 5 ' 2 near £ = and as a result, 
the density behaves as < I )4 oc r~ 4 . The self-consistent r c models and the J m models (see the 
following section) presented here are infinite-radius and finite-mass models with an r -4 density 
profile. We have checked that the finite and non-zero N{£ = 0) property also holds for the angular 
momentum restricted density of states. 



5. /o = l/T r model 

The numerical solution of Poisson's equation for /o = 1/T r follows closely the analytics in §pl 
The areas of integration in (£, J 2 ) space are determined by ( |i~2| , |l8| ). Starting with the initial guess 
of the potential given in appendix |A| or corresponding one for angular momentum cutoff, the radial 
period T r (£,ro) is calculated by a root solving and integration routine within the bound given by 
£ < £f(ro,ro) for tq < r c and £ < £f(r c ,ro) for vq > r c in case of the r c model. Similarly, the 
bound £ < £f(ro,ro);ro < r max (£) was used in the case of the J m model where r max (£) is radius 
of largest allowed circular orbit. A lookup (interpolation) table for T r (£,ro) was prepared. For 
a given radius, r, the regions of integration are determined by A or B and then the distribution 
function 1/T r is integrated to calculate the density profile using 
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1 d , 2 d$, 
(r — ) = 


- A 


X(r;T r 1 ) r c model 


1 


K(r;T~ l ) J m model 



(26) 



T r (S,r) = 2 / dr (2(-£-$)-J 2 /r 2 )- ' , (27) 

Jr p 

where r p and r a are the turning points and the value of K is chosen so that the total mass is 1 
and G = 1. The next iterate of the potential is then trivially obtained from 

/•oo 

H0 = -J C 2 M(0dH (28) 

where «M(£) is the mass fraction inside £ that is calculated from the density, the RHS of (p6|). 
Here, $(oo) is taken to be zero to be consistent with the initial guess of £~ 3 > 2 for the radial period. 

The numerical code was used to verify the £ 3 ' 2 solution (the initial guess) and vice versa. 
A high precision routine was used to calculate the radial period for an arbitrary potential and 
this was tested against the well-known forms for isochrone and Kepler potentials and found to be 
precise to within a millionth of the correct value. The verification of the analytics established the 
robustness of the numerics used. The numerical scheme converges rapidly to a solution in a few 
iterations. This indicates that the properties of these models are close to those of the £ 3 ' 2 models. 
Fig. [7] illustrates that the r c model is nearly isochronic. For many purposes, one could have been 
satisfied with the £ 3 ' 2 models, which are much easier to implement than the 1/T r models, for 
exploratory work. But this is strictly with the hindsight provided by our constructing the models 
in the first place. 



5.1. The energy distribution 

Let us now look at the energy distribution, N(£). For a spherical system with the distribution, 
f(£ , J 2 ), one can write the differential energy distribution (see Binney &: Tremaine, eqn 4P-11) as 

N(£) = 4ir 2 f f(£, J 2 )T r (£, J 2 ) dJ 2 . (29) 

In the case of / = A C(rj)/T r , we obtain 

N(£) = 4ir 2 A fc(ri)dJ 2 = Air 2 A J 2 {£) (30) 

where J~(£) is the constraint boundary in the £ — J 2 plane. In the case of the r c model, the shape 
of the constraint in this plane (c.f. Kl, § |3.1| ) determines the energy distribution which is given by 



r 2 A j t 3 $'(t); where -$(£) -W(t)/2 = £ £ > £ c 
-2(£ + $ c )r c 2 £ <£ c 



N c {£)=A<k 2 A{ v Zo.~?2 ' "V"*- 1 - Y/: (31) 
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where £ c is the energy of a circular orbit at r c as defined in §0, and similarly for J m model (c.f 



3.2), we obtain 



2A jt 3 $'(t); where - $(t) - t$'(t)/2 = £ £ > £ n 

J 2 £ < £ 



N j (£) = 4n 2 A\ " 2 , w ' «"— W--W— ---m (32) 



where the energy £ m of the largest allowed circular orbit is given by 

£ m = -*(*) - iV(i)/2; £ 3 $'(£) = «7^. (33) 

The Keplerian limit of eqn ( |3lD and eqn p2] ) agrees exactly with the Keplerian limits of £ 3 ' 2 
counterparts, when one takes into account the fact that for a Kepler potential, T r = ttGM£~ 3 ' 2 /v2. 
Therefore A/A' = 7r/\/2 in units where GM = 1. It is clear that, near small £ , the form of N(£) 
is linear for the r c model and the slope is —8ir 2 Ar 2 (positive in N(E)) while the N(£) is flat for 
the J m model. From the dynamical arguments in §g we expect orbits with small T r to be more 
populated than the ones with larger radial periods. However, we have no reason to adhere to the 
strict 1/T r form for small T r . It is worth checking how sensitive our results are to a modification 
in which /q is flatter at small T r while retaining the 1/T r asymptotic behavior. We thus try 

f Q (E,J 2 ) = -^— (34) 

i r + J- 

where To is a parameter which is chosen to be of the order of the radial period of the harmonic 
oscillator, T^, near the bottom of the well and given by 

T ^wm (35) 

Again, the maximum deviation from the initial guess upto the second iterate in the density was 
found to be less than 1% for To = 0.5, 1,2 T^. This indicates that the distribution function is 
probably stable and insensitive to small changes in /o but sensitive to the choice of the constraint, 
C(rj). This and the similarity of the £ 3 ' 2 models to the 1/T r counterparts can be explained by the 
fact that /o is dominated by £ 3 / 2 (see the contour plot of T r £ 3 / 2 , Fig ^). 



6. Properties of the solutions and comparison with simulations 

We discuss the analytic results in the preceding sections and the properties of the r c and 
J m solutions and compare it with relevant simulations and the r 1 ' 4 law. As mentioned in §ffl the 
simulations of direct relevance are cold collapse simulations, in particular the C runs of van Albada 
(1982) who has presented the density, N(£), velocity dispersion and anisotropy profiles of the final 
configurations. In our figures |l], H|, ||, || and g a Henon isochrone potential was used to illustrate 
the allowed region in phase space. 

Scaled quantities and their physical units 
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The models have two free parameters, the total mass M, and either the cutoff radius r c or a 
maximum value of angular momentum, J m . Below, we compare the r c models with the J m models. 
The scale radius, r s is the scale radius of the J m model and equals the core radius, r c , in the r c 
model. If the total mass and scale radius are set to unity, the model has no free parameters. 

The potential in the r c solution is scaled according to 

GM GM , % 

$ = 6 = 0.265 9 (36) 

a r c r c 

and similarly the scale for potential in the J m model is 0.462 GM/r s . The density scales as 
(l/47ra) Mjr\ = 0.021 Mjr\ for the r c model whereas it scales as 0.0368 Mjr\ for the J m model. 
The radial orbit period can be written as 

T - *° I' 2 ** (37) 



and hence we define a unit, T c = 2r c /y / a = 2y/a \Jr\jGM. For the r c model, a = 3.38 (T c = 3.68) 
and a = 2.16 (T c = 2.94). 

The density profile 

For large r, the densities, p(r), for both the models scale as r -4 ; this has been derived 
analytically in appendix [A| and is also apparent from Fig. [81 The density has a sharp break at 
r = r c for the r c models. The fraction of mass in the core, is about 5% for both models. The 
density continues to decrease gently with r upto about 5r c , beyond which the slope changes rather 
abruptly to a much steeper value, and rapidly converges to the asymptotic r -4 profile (see the 
log- log plot in the lower panel of Fig. ph. The radius at which the break occurs is a few r s (see eqn 



( A14 )) for both models. This is seen in the C runs of van Albada (1982), see Fig 6 in the paper. 
A similar break is also seen in the cosmological simulations of Tormen, Bouchet & White (1997) 
and Moore et al. (1998) where the power law is roughly r _1 for the inner region and between r~ 3 
and r -4 for the halo. 

The differential energy distribution 

Here we use E = —E in order to compare our results with simulations. The differential 
energy distribution, N(E), is a useful indicator of the phase space structure. The differences 
in the analytic forms of N(E) for the different constraints are indicated by eqns (p5|, 31, & p2]) 



which corroborate each other in the Keplerian limit (this is a reflection of the shape of the C(r/) 
in the £ — J 2 plane). Whereas the real space structure (especially at large r) for the two models 
are similar, the plots in Fig. |9| demonstrate how different the models are at high energies; N(E) 
for the J m model rises sharply, reaches a maximum at some E, and thereafter is independent of 
E. For the r c model, N(E) is smooth, rising gradually, and becomes a strictly linear function 
of E (see Fig 7 of van Albada 1982, runs C2 and C3 for a linear plot). The linear region in our 
model begins at energies 0.7 of the well depth (c.f appendix^), i.e, it is dominantly linear and 
is consistent with C2 and C3. The log plot, Fig. is also consistent with Fig 4-20 in Binney &; 
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Tremaine of run C3 and the simulations of Spergel & Hernquist (1992), their Fig. 2. Also it is 
clear from eqn (^) that the non-zero intercept of the constraint, C(rj), which is J~(£ = 0), has 
the desirable effect of implementing Jaffe's (1987) insight that the differential energy distribution 
N(E) tends to a non-zero value as E — ► from below, since the probability of ejection from the 
core is not expected to be sensitive to small changes in the final energy. 

De Vaucoleurs r 1 ' 4 law 

In both cases, the surface densities provide reasonably good fits to the de Vaucoleurs r 1 ' 4 law 
in the range 0.1-8 r e as indicated in Fig. [To], assuming a constant mass-to- light ratio. This includes 
the range in radii (0.1 -2 r e ) over which the r 1 ' 4 law provides excellent fits to the brightness 
profiles of ellipticals (Burkert 1993). The slope in the figure is close to —8 (the standard slope is 
-7.67) where £ was normalized to its value at the center unlike the usual normalization by the 
extrapolated peak value. 

The anisotropy profile 

The DF, eqn (|J), naturally gives radially anisotropic models where most of the mass is outside 
the core radius. The anisotropy profile and velocity dispersions of the £ 3 ' 2 models which is very 



similar to that of the T r models. Fig. 11 shows a plot of the run of the anisotropy parameter, 
(3, with radius. Both the r c and J m models are very nearly isotropic within the core and rapidly 
become anisotropic. The Fig. ^ shows the radial velocity dispersions which indicate that the r c 
model is slightly more anisotropic than the J m model. The dispersion and anisotropy profiles are 
very similar to the Fig 8 of van Albada (1982). 



7. Discussion 

In §[y we made an assumption, based on a plausible picture of diffusion in action space, 
that A(E) = constant; as discussed in the previous section, the agreement with simulations is 
encouraging. However, the physical picture of mixing without friction leads to a constant A only 
if the constraints are hard and the mixing lasts enough time, two things we cannot lay down a 
priori. So the success of the simplest such model does teach us that this picture is a good first 
approximation, possibly improvable. 

We now compare and contrast our approach with two earlier papers which can be regarded 
as fitting into the same broad framework of constrained violent relaxation. MTJ explored DFs 
with Gaussian and Lorentzian cutoffs in J, multiplied by an exponential function of E; thus 
their models have one more parameter, a temperature, 1/(3 (c.f. Stiavelli & Bertin 1987). The 
infinite temperature limits of their models, in common with our sharply cutoff J m model, have 
an asymptotically flat N(E) near E = 0. In order to obtain an increasing N(E) with the MTJ 
models, one requires a negative (3. As discussed earlier, our r c model gives an increasing N(E) 
without the need for an exponential factor. We therefore find the r c model is preferable to the J m 
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model. We note that even if the cutoff function, C(rj) is not as abrupt as given in eqn (||), N(E) is 
expected to behave in a similar fashion. 

SH proposed a kinetic model of the wave-particle interaction process, with an associated 
variational principle analogous to Boltzmann's H-theorem, and made a comparison with their own 
extensive simulations. Their description of energy changes occurring by a sequence of kicks is in 
fact close to the picture here. We used this to motivate a constant phase space density modulated 
by a pericenter or angular momentum cutoff and Tremaine's incompleteness factor (T e /T r ). 
In contrast, SH assume that the kinetics drives the system to the maximum of their entropy 
functional. However, their Boltzmann-like factor contains a negative temperature (as in the MTJ 
models), but this equilibrium is meaningful only in systems for which the density of accessible 
states decreases rapidly with energy (c.f. § 73 of Landau & Lifshitz 1980). We suggest that the 
SH kinetic picture, with which we are in broad agreement, would not actually lead to a negative 
temperature distribution. The more appropriate physical picture is one of mixing in phase space. 



8. Conclusions and Caveats 

We have presented a semi-analytic model of violent relaxation that includes a new picture 
of diffusion in phase space with a novel implementation of the pericenter constraint, and the 
1/T r incompleteness factor. Notions of even a partial thermal equilibrium with Boltzmann-like 
exponential factors play no role - a property we regard as a virtue in describing a collisionless 
system. The rise in the energy distribution function which such factors (with negative 
temperatures!) mimicked in earlier work now arises naturally from the pericenter cutoff in our 
calculation. The resulting properties of density, surface brightness, energy distribution, anisotropy 
profiles are in good agreement with simulations. The fact that the properties of the models are 
parameter free (M and r s are merely used to normalize), may be a considered a virtue, as they 
demonstrate that the constraints explored here seem to capture the essential details in the case 
of cold non-rotating collapse. More realistic systems can be obtained if one deviated from the 
simplifying assumptions and comprehensively explored associated models such as l/(T r + To) in 
§||, albeit with more parameters. Even though sophisticated numerical codes that now exist to 
perform N-body experiments, semianalytic models help in understanding their output. They also 
have the advantage of fitting simulations and observed galaxies reasonably well. Possible directions 
of future work include a more detailed comparison of our models with numerical simulations, 
investigations of stability and extension to axisymmetric systems. 



A. £ 3//2 model with pericenter cutoff 

We have carried out a detailed study of models with /o, the first factor in eqn (y), chosen 
to be £ 3 ' 2 , with both kinds of cutoff, viz. pericenter and angular momentum. Exisiting work, 
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by MTJ, uses a gaussian rather than sharp cutoff in angular momentum. The pericenter cutoff 
has not been implemented previously and §^ gives some details to enable the interested reader 
to see how the sharply cutoff models are constructed, and other details of its properties along 
with comparisons to simulations can be found elsewhere (Mangalam & Sellwood 1999). As a first 
approximation and the simplest model, we consider the distribution function /o = £ 3 ' 2 , with a 



sharp pericenter cutoff at r c . Then the integral in §3.1 can be written as 



h(r 



1 rJl fB{r,J 2 ) 

£3/ 2) = _L_j e dJ 2J V &£ £*/\-£ + BY 1 ' 2 

= ^|(-*) 8 (Al) 



Similarly, the integral, 

but where J 2 is the lower limit to J 2 and 



l_(r;£ 3 / 2 )^l_ = —£ 3 (A2) 



/C-(r;^/ 2) = ^ ( _„_^)3 (A3 ) 



where J^ is the lower limit to J 2 . Now the "+" integral works out to be 

J + (r;5 3 / 2 ) = -4/ d££ 3 /\(-<I>-£)2ri-j2 

r l Jo v J j2 

= JL£3 (1 _^ (r 2_ r 2 )/r 2) (A4) 



Gathering all the integrals in X, defined above and (12), and absorbing the numerical factor into a 



positive constant K, Poisson's equation can be written as 



1 d 'V 2x__ rr) -® r ^ r c 



($> r *) = K\ (r2 _ r2) _ 5 /2 (A5) 

This model is a polytrope of index 3 for r < r c and we recover the standard results for a 
distribution function of the form £ 3 ' 2 without the pericenter cutoff (See Binney & Tremaine 1987, 
eqn (4-108c)). For r > r c , the contribution to the density is dominated more by orbits with largely 
radial velocities and the first of the 2 terms limits the orbits to those which have angular momenta 
less than the bound specified by Ai- In the limit of r c = 0, the RHS of ( |A5[ ) for r > r c vanishes as 
expected. We now scale the radius by A = r/£, and the potential by a = &/6, similar to scalings 
employed for the polytrope problem (with the difference that we have the normalization constant 
K) and obtain the following relations 

Ka 2 X 2 = 1 (A6) 

GM/(\a)=a = -J o fe 3 d£ - J £ {-(£ 2 - £ 2 )- 5/ W - 0£ c f + ^ 3 } d£ (A7) 
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where a is a geometric factor that depends on the solution. Here we pick X/r c = £ c = 1, (so that 
a = -^-y and K = (^jy) 2 ) such that the cutoff radius in these units is unity and it simplifies the 
task of finding the solution. This model has two parameters, the total mass M and the cutoff 
radius r c and a unique solution can be found from the resulting equation, where 6 C = 6(1) = 6\ 

i d /Q/t2 , f -o 3 e<i 



e dc {H) -\ ^f^{ee - e x f - & 3 z > i (A8) 



with the boundary conditions 



?'(0) = & 0(oo) = 0. (A9) 



The former is a consequence of an implicit assumption of a non-existence of a central point 
mass and the latter boundary condition is enforced to be consistent with the distribution 
/o = £ 3 ' 2 ~ 1/T r vanishing at £ = 0. It is interesting to study the asymptotics of this model by 
changing to a convenient variable u = l/£. For u < 1, we obtain 

u 4 9"(u) = {9- 9 lU 2 f{l - u 2 )- 5 / 2 - 9 3 (A10) 

As r — > oo or near u = 0, we can expand in powers of u to obtain 

u i 0"(u)/9 2 = -3u 2 9 2 + 3£u 4 + (5/2)£V - u e + 0(u 7 ) (All) 

where = 6/6\, which leads to the following solution of 9(u) near u = 

d(u) = cm - (3/2)cl0 2 u 2 + 0(u 3 ) (A12) 



where c\ = a/^i since $ — ► —GM/r c u or — > —a u. As a result the density, p, in units of M/r 3 
asymptotically behaves as 

p = u A 9"{u) = ^{-3c 2 u 4 + [3ci-9c 2 fl 2 + (5/2)c?]u 5 } + 0(u 6 ). (A13) 

Hence, as r — > oo, p ex r~ 4 . Now the radius beyond which the leading term dominates, which we 
call as the break radius, is defined as 



n 



- + U- w 2 

C\ D 



(A14) 



Next, we solve the equations, (|A8|) and (A9). It is convenient to use Lane-Emden solutions 



in the region, (0 < £ < 1), and use ( |A10[ ) in the region 1 > u > 0, with the boundary condition 
9(u = 0) = 0. A sufficiently accurate solution for the region < £ < 1 can be obtained by power 
series expansions. The solution upto order 12 is given in terms of the well depth 6(0) = 6q can be 
obtained from power series expansions. Using that solution for a trial well depth, 6q, 6(u = 1) and 
6'(u = 1) = — #'(£ = 1) are calculated and ( |A10| ) is numerically integrated to 6(u = 0). A value 
of 6q is for which the boundary condition, 6(u = 0) = is satisfied is determined iteratively. This 
occurs for 

6 = -0.821 6'(u = 1) = -0.15264 6(u = 1) = -0.737234 (A15) 

and the corresponding 9(£) is shown in Fig. ^. Note that for this solution, a = 3.38. 
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Fig. 1. — The two possible locations of r c , r c i < tq and r c i > vq are indicated in the figure where 
vq is the circular orbit radius. The allowed orbits are those for which £ < £f(r c ,ro), as indicated 
in the figure, and £ < £/(ro,ro), the energy of the circular orbit, respectively. 

Fig. 2. — The applicable bound for r < r c is the circular orbit energy, £f(ro,ro) and for r > r c , it 
given by £ < £f(r c , ro) as explained in Fig. 1. Lines of turning points are drawn for r > r c , r < r c 
and r = r c . 

Fig. 3. — The areas of integration enclosed by the bold lines-the triangle shaped Ai, for r < r c and 
the wedge shaped A2, for r > r c are indicated in the figure. The limits used in evaluating the 2 
integrals, £*, J 2 and J m are also indicated. The £ intercept is — <&(r) and the foot of A2 is — 2<J> c r 2 . 

Fig. 4. — The solution for the potential in the £ 3 ' 2 — r c model is shown in the figure in units of 
0.265 GM/r c where r c is the unit of radius. 

Fig. 5. — Jmax(£> r ) f° r the case £ > £ c and r > r c is at the point of intersection. In this case, 
£ > £*, since £ c > £*. This implies that J^ax = ^e = ~~ 2(£ + <l>)r 2 . When r < r c , the pericenter 
constraint is satisfied and hence J^ax = Je- 

Fig. 6. — J;^ aa ,(<? > £ c , r) for £* > £ and £* < £ are shown in the figure. £ c > £ > £*% for a radius 
r = X2 > t c and J^aa; = — 2(f + <&)r 2 , given by the point of intersection at the energy £ > £ c . 
Similarly, £ < £*i < £ c for r = x\ > r c and J^ nax = — 2(£ + ^c)^- Also r c < x\ < rj < X2 where 
ri, given by {£ — &i)r] = {£ — $ c )^c) i s the apocenter of an orbit with energy — £ and pericenter 
r c - 

Fig. 7. — The top panel shows a contour plot of the variation of T r £ 3 ' 2 in the allowed £ — J 2 plane. 
The flat contours indicate that T r is very nearly isochronic and close to £~ 3 ' 2 . The deviation is 
depicted in the lower panel which shows a section taken at J 2 / ' J 2 ir = 0.5, where J 2 ir (£) is the 
angular momentum of a circular orbit at a given energy, £. 

Fig. 8. — The upper panel shows the density as a function r/r s , where r s = r c for the r c models, 
and r s is the scale radius of the J m models. The former shows a sharp drop at r c due to the 
pericenter constraint whereas the latter has a smooth profile. The log-log plot in the lower panel 
shows a break near 5 r c and an asymptotic form of r~ 4 from 10 -100 r c for both models. The 
density is shown in units of 0.021 Mjr\ for the r c model and 0.0368 Mjr\ for the J m model. 

Fig. 9. — The N(E) for the J m model rises abruptly and then is flat, whereas the N(E) for the r c 
model increases gradually. The well depth in both cases was taken to be 1. 

Fig. 10. — The surface density, S is an excellent fit to the r 1 ' 4 law for 0.1 < r/r e < 8. r e = 6.65r c 
for the r c models (lower curve), and r e = 4.03r s for the J m models (upper curve). 

Fig. 11. — A plot of the anisotropy parameter, (3(r) = 1 — v 2 /(2 v 2 ), for both models. The core is 
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nearly isotropic and becomes nearly radially anisotropic at 10 r s . The r c model (upper curve) is 
slightly more anisotropic than the J m model (lower curve). 

Fig. 12. — A plot of the normalized radial velocity dispersion, a^(r) = v^(r)/v^(0), for both 
models. The r c model (lower curve) is slightly more radially anisotropic than the J m model (upper 
curve) . 
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